単日データを処理するための関数を定義しておきます。
| 関数名 | 処理概要 | 必要なパッケージ |
|---|---|---|
lagdiff(n) |
単日データの前日差を求める | dplyr |
ma7(n) |
単日データの7日間(1週間)移動平均を求める | zoo |
ma28(n) |
単日データの28日間(4週間)移動平均を求める | zoo |
ma(n, k) |
単日データのk日間移動平均を求める | zoo |
lagdiff <- function(n) { # 前日差を求める関数
n - dplyr::lag(n, default = 0L)
}
ma7 <- function(n) { # 移動平均(7日)を求める関数
zoo::rollmeanr(n, k = 7L, na.pad = TRUE)
}
ma28 <- function(n) { # 移動平均(28日)を求める関数
zoo::rollmeanr(n, k = 28L, na.pad = TRUE)
}
ma <- function(n, k = 7L) { # 移動平均(k日)を求める関数
zoo::rollmeanr(n, k = k, na.pad = TRUE)
}
累計は Base R の cumsum 関数を利用します。
その他、グラフで繰り返し使う文字列を変数として定義しておきます。
subtitle <- paste0("Generated @", lubridate::now())
caption <- "Data Source: covid19japan.com"
個票データの集計に限らず、データをインポートした際には目的に応じて各変量(変数)の変数型を変更します。特に水準ごとに層別処理を行いたい場合には因子型に変換しておくと便利です。また、結合したいデータと名称や体系を合わせておくこともポイントです。
都道府県に関するデータは Gist で公開しているデータを用います。データ都道府県コード順に並んでいますので、この並びを利用して forcats パッケージを用いて因子化します。Base R の as.factor 関数を用いると水準がアルファニューメリック順に並べ変えられてしまう点に注意してください。
prefs <- "https://gist.githubusercontent.com/k-metrics/9f3fc18e042850ff24ad9676ac34764b/raw/f4ea87f429e1ca28627feff94b67c8b2432aee59/pref_utf8.csv" %>%
readr::read_csv() %>%
dplyr::mutate(
# Googleの予測データと結合するためにコード体系を合わせる
japan_prefecture_code = paste0("JP-", `コード`)
) %>%
dplyr::select(
# Googleの予測データと結合するために名称を変更する
japan_prefecture_code, prefecture_name = pref,
# 日本語の変数名は扱いにくいので英語名に変更する
pref = `都道府県`, region = `八地方区分`, pops = `推計人口`
) %>%
dplyr::mutate(
# 水準ごとに表示させるために因子化する(あらかじめデータをコード順に
# 並べておくことが因子化の際のポイントのひとつ)
japan_prefecture_code = forcats::fct_inorder(japan_prefecture_code),
pref = forcats::fct_inorder(pref),
region = forcats::fct_inorder(region),
pops = as.integer(pops)
)
prefs
個票データは Covid19Japan.com のデータリポジトリ で公開されているデータを用います。データは参照しているソースが更新されると自動的に更新されますので、取得タイミングによっては全てのデータが揃っていない可能性があります。
JSON形式データの各項目については データフォーマット を参照してください。
df <- "https://raw.githubusercontent.com/reustle/covid19japan-data/master/docs/patient_data/latest.json" %>%
# JSON形式をデータフレームとして読み込む
jsonlite::fromJSON() %>%
dplyr::select(patientId, date = dateAnnounced, gender,
pref = detectedPrefecture, patientStatus, knownCluster,
confirmedPatient, ageBracket,
deceasedDate, deceasedReportedDate) %>%
dplyr::filter(confirmedPatient == TRUE) %>%
dplyr::mutate(
date = lubridate::as_date(date),
gender = forcats::as_factor(gender),
pref = stringr::str_to_lower(pref),
patientStatus = forcats::as_factor(patientStatus),
cluster = dplyr::if_else(!is.na(knownCluster), TRUE, FALSE),
ageBracket = forcats::as_factor(ageBracket),
deceasedDate = lubridate::as_date(deceasedDate),
deceasedReportedDate = lubridate::as_date(deceasedReportedDate)
) %>%
# 都道府県データと結合
dplyr::left_join(prefs, by = c("pref" = "prefecture_name")) %>%
dplyr::select(-pref) %>%
dplyr::rename(pref = pref.y) %>%
# 因子型の欠損値を水準化する
dplyr::mutate(
japan_prefecture_code = forcats::fct_explicit_na(japan_prefecture_code,
na_level = "JP-48"),
pref = forcats::fct_explicit_na(pref, na_level = "空港検疫"),
region = forcats::fct_explicit_na(region, na_level = "空港検疫"),
gender = forcats::fct_explicit_na(gender, na_level = "非公表"),
ageBracket = forcats::fct_explicit_na(ageBracket, na_level = "非公表"),
patientStatus = forcats::fct_explicit_na(patientStatus,
na_level = "Unknown")
)
df
作成したデータフレームを確認します。意外と欠損値(n_missing項参照)が多いデータであることが分かります。
df %>%
skimr::skim()
| Name | Piped data |
| Number of rows | 429860 |
| Number of columns | 14 |
| _______________________ | |
| Column type frequency: | |
| character | 2 |
| Date | 3 |
| factor | 6 |
| logical | 2 |
| numeric | 1 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| patientId | 0 | 1.00 | 1 | 16 | 0 | 429860 | 0 |
| knownCluster | 427355 | 0.01 | 3 | 88 | 0 | 233 | 0 |
Variable type: Date
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| date | 0 | 1 | 2020-01-15 | 2021-02-26 | 2020-12-26 | 395 |
| deceasedDate | 429480 | 0 | 2020-02-13 | 2020-11-19 | 2020-05-08 | 151 |
| deceasedReportedDate | 429530 | 0 | 2020-02-13 | 2020-10-17 | 2020-05-16 | 131 |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| gender | 0 | 1 | FALSE | 3 | 非公表: 320188, M: 61370, F: 48302 |
| patientStatus | 0 | 1 | FALSE | 9 | Unk: 427326, Hos: 1261, Dec: 372, Hom: 315 |
| ageBracket | 0 | 1 | FALSE | 13 | 非公表: 320285, 20: 29433, 30: 19042, 40: 16089 |
| japan_prefecture_code | 0 | 1 | FALSE | 48 | JP-: 110240, JP-: 46965, JP-: 44660, JP-: 29174 |
| pref | 0 | 1 | FALSE | 48 | 東京都: 110240, 大阪府: 46965, 神奈川: 44660, 埼玉県: 29174 |
| region | 0 | 1 | FALSE | 9 | 関東地: 224459, 近畿地: 83404, 中部地: 43186, 九州地: 37312 |
Variable type: logical
| skim_variable | n_missing | complete_rate | mean | count |
|---|---|---|---|---|
| confirmedPatient | 0 | 1 | 1.00 | TRU: 429860 |
| cluster | 0 | 1 | 0.01 | FAL: 427355, TRU: 2505 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| pops | 2232 | 0.99 | 7917.71 | 4229.68 | 560 | 5107 | 7537 | 13822 | 13822 | ▆▅▆▆▇ |
最初に単純集計を行ってみます。
全国集計は日付データ(data)に基づく日次の単純集計です。厚生労働省のオープンデータにおける陽性者数データに相当します。単純集計は dplyr::group_by と dplyr::summrise を用いることで簡単に計数できます。ただし、データがない日は集計結果として現れませんので tidyr::complete で補完しておく必要があります。
その後、集計結果を元に前日差、累計、移動平均を計算します。
japan_daily <- df %>%
dplyr::group_by(date) %>%
dplyr::summarise(n = dplyr::n()) %>%
dplyr::ungroup() %>%
tidyr::complete(
date = seq.Date(from = min(date), to = max(date), by = "day"),
fill = list(n = 0L)
) %>%
dplyr::mutate(
diff = lagdiff(n), # 前日差
cum = cumsum(n), # 累計
ma7 = ma7(n), # 移動平均(7日)
ma28 = ma28(n) # 移動平均(28日)
)
japan_daily
| 変数 | 変数の意味 | 変数型 | 備考 |
|---|---|---|---|
| date | 陽性判定者の報告日 | Date | 陽性判定日とは異なる |
| n | 陽性判定者数 | Integer | |
| diff | 同前日差 | Integer | 初日は陽性判定者数に同じ |
| cum | 同累計 | Integer | |
| ma7 | 同移動平均(7日) | Double | |
| ma28 | 同移動平均(28日) | Double | 4週間(1ヶ月相当) |
dplyr::group_by で指定する変数を変えることで各変数の水準ごとの集計ができます。これは、前述の skimr::skim によるサマリーと同値になりますので下記以外の集計は省略します。
df %>%
dplyr::group_by(gender) %>%
dplyr::summarise(n = dplyr::n()) %>%
dplyr::ungroup()
クロス集計は日付と地方区分、日付と都道府県などのように複数の変数における水準ごとの集計です。tableなどを利用しても集計可能ですが、ggplot2 パッケージによる可視化を考慮し、ここではデータフレームを用いた集計を行います。
以下のような構造を取る個票データをクロス集計するとします。
df <- data.frame(
date, # Date型の日付データ
key, # グルーピングのためのデータ(例:都道府県、地方区分など)
... # 順不同
)
基本的には単純集計と同じで dplyr::group_by と dplyr::summarize で集計します。補完も単純集計と同じく tidyr::complete で行うことができます。
x <- df %>%
dplyr::group_by(date, key) %>% # クロス集計対象を指定する
dplyr::summarise(n = dplyr::n()) %>% # n()は個数をカウントする関数
dplyr::ungroup() %>% # 最後にungroupするのがポイント
tidyr::complete( # 暗黙の欠損を補完する
date = seq.Date(from = min(date), to = max(date), by = "day"),
key, fill = list(n = 0L) # 個票がない=陽性者ゼロ
)
クロス集計後に前日差や累計などを算出するには tidyr::nest,tidyr::unnest と purrr パッケージを組み合わせた処理が確実かつ高速です。
x %>%
dplyr::group_by(key) %>% # 上記の日次集計を元に各種計算を行う場合
tidyr::nest() %>% # nest & unnest 処理を利用すると確実
dplyr::mutate( # ネストしたデータはpurrrで処理する
diff = purrr::map(data, ~ lagdiff(.$n)), # 前日差
cum = purrr::map(data, ~ cumsum(.$n)), # 累計
ma7 = purrr::map(data, ~ ma7(.$n)), # 移動平均(7日)
ma28 = purrr::map(data, ~ ma28(.$n)) # 移動平均(28日)
) %>%
tidyr::unnest(cols = c(data, diff, cum, ma7, ma28))
purrr パッケージの使い方がよく分からない場合は以下の方法でも可能です。
x %>%
dplyr::group_by(key) %>% # ここがポイント
dplyr::mutate(
diff = lagdiff(n),
cum = cumsum(n),
ma7 = ma7(n),
ma28 = ma28(n)
) %>%
dplyr::ungroup()
前述のクロス集計処理を daily_aggregate 関数として定義しているので、これを用いて集計を行っています。
region_daily <- df %>%
daily_aggregate(date, region) # 前述の集計処理を関数化したもの
region_daily
地方区分集計と同様に集計します。
pref_daily <- df %>%
daily_aggregate(date, pref)
pref_daily
ageBracket_daily <- df %>%
daily_aggregate(date, ageBracket)
ageBracket_daily
年代別の地方区分別集計ならびに都道府県別集計は割愛します。理由はVisualize章のグラフを見て頂ければ分かると思います。
cluster_daily <- df %>%
daily_aggregate(date, cluster)
cluster_daily
クラスター別の地方区分別集計ならびに都道府県別集計は割愛します。理由はVisualize章のグラフを見て頂ければ分かると思います。
作成した集計データを可視化します。
japan_daily %>%
dplyr::filter(date == max(date))
単日と累計が二桁違うので、縦二軸のグラフを描きます。
title <- "【全国】陽性者数(単日)"
xlab <- ""
ylab <- ""
sec_scale <- 50 # 縦二軸のスケーリング
japan_daily %>%
ggplot2::ggplot(ggplot2::aes(x = date)) +
ggplot2::geom_bar(ggplot2::aes(y = n), stat = "identity", width = 1.0,
fill = "dark gray", alpha = 1.0) +
ggplot2::geom_line(ggplot2::aes(y = ma7), linetype = "solid",
colour = "gray10", size = 0.5) +
# 第二軸を利用するグラフを描画する際はスケーリング調整する必要あり
ggplot2::geom_line(ggplot2::aes(y = cum / sec_scale),
colour = "dark green", size = 1.0) +
# 横軸表示の指定
ggplot2::scale_x_date(date_breaks = "1 month", date_labels = "%y/%m") +
# 二軸表示のための軸属性の指定
ggplot2::scale_y_continuous(
# 第一軸のラベル(スケールは自動調整)
name = "陽性者数(灰)・移動平均(黒)",
# 第二軸の指定(第一軸にあわせたスケーリング)
sec.axis = ggplot2::sec_axis(~ . * sec_scale,
name = "累積陽性者数(濃緑)")) +
# 縦軸の装飾
ggplot2::theme(
axis.text.y.left = ggplot2::element_text(colour = "gray10"),
axis.line.y.left = ggplot2::element_line(colour = "gray10"),
axis.text.y.right = ggplot2::element_text(colour = "dark green"),
axis.line.y.right = ggplot2::element_line(colour = "dark green")
) +
# グラフ全体のラベル指定
ggplot2::labs(title = title, subtitle = subtitle, caption = caption,
x = xlab, y = ylab)
以降、描画コードは割愛しますが、このように個票を任意の変数の水準ごとに集計すると様々なデータの見かたができるようになります。
japan_daily %>%
dplyr::mutate(weekday = lubridate::wday(date, label = TRUE),
weeks = lubridate::epiweek(date)) %>%
dplyr::mutate(weeks = dplyr::if_else(date > "2021-01-02",
weeks + max(weeks), weeks)) %>%
# print() %>%
ggplot2::ggplot(ggplot2::aes(x = weekday, y = weeks, fill = n)) +
ggplot2::geom_tile() +
ggplot2::scale_fill_gradient(low = "white", high = "red") +
# ggplot2::coord_equal() +
# ggplot2::coord_flip() +
ggplot2::labs(title = "陽性者数ヒートマップ", x = "", y = "",
subtitle = subtitle, caption = caption)
japan_daily %>%
dplyr::filter(n > 0) %>%
dplyr::mutate(weekday = lubridate::wday(date, label = TRUE)) %>%
ggplot2::ggplot(ggplot2::aes(x = weekday, y = n)) +
ggplot2::geom_violin(ggplot2::aes(fill = weekday), alpha = 0.25) +
ggplot2::stat_boxplot(geom = "errorbar", width = 0.1) +
ggplot2::stat_summary(ggplot2::aes(group = 1, colour = weekday),
fun.y = median, geom = "line") +
ggplot2::geom_boxplot(width = 0.1) +
ggplot2::labs(title = "曜日別陽性者数分布", x = "", y = "",
subtitle = subtitle, caption = caption)
japan_daily %>%
dplyr::mutate(month = dplyr::if_else(lubridate::year(date) > 2020,
lubridate::month(date) + 12,
lubridate::month(date))) %>%
dplyr::filter(n > 0) %>%
dplyr::mutate(weekday = lubridate::wday(date, label = TRUE)) %>%
ggplot2::ggplot(ggplot2::aes(x = weekday, y = n)) +
ggplot2::geom_violin(ggplot2::aes(fill = weekday), alpha = 0.25) +
ggplot2::stat_summary(ggplot2::aes(group = 1, colour = weekday),
fun.y = median, geom = "line") +
ggplot2::geom_boxplot(width = 0.1) +
ggplot2::facet_wrap(~ month, nrow = 12, scales = "fixed", dir = "v") +
ggplot2::labs(title = "曜日別・月別陽性者数分布", x = "", y = "",
subtitle = subtitle, caption = caption)
subset <- region_daily %>% dplyr::mutate(key = region)
title <- "【地方別】陽性者数分布"
xlab <- ""
ylab <- ""
ncol <- 3
subset %>%
dplyr::filter(n > 0) %>%
ggplot2::ggplot(ggplot2::aes(x = reorder(key, n, median),
y = n, colour = key)) +
ggplot2::stat_boxplot(geom = "errorbar", width = 0.3) +
ggplot2::geom_boxplot() +
ggplot2::theme(legend.position = 'none') +
ggplot2::coord_flip() +
ggplot2::labs(title = title, x = "", y = "",
subtitle = subtitle, caption = caption)